Effective Magnetic Fields in Graphene Superlattices 



Jianmin Sun and H.A. Fertig 
Department of Physics, Indiana University, Bloomington, IN 47405 



o 

(N 

C 
3 



L. Brey 

Institute) de Ciencia de Materiales de Madrid ( CSIC), 

(Dated: June 8, 2010) 



Cantoblanco 28049, Spain 



We demonstrate that the electronic spectrum of graphene in a one-dimensional periodic potential 
will develop a Landau level spectrum when the potential magnitude varies slowly in space. The 
effect is related to extra Dirac points generated by the potential whose positions are sensitive to its 
magnitude. We develop an effective theory that exploits a chiral symmetry in the Dirac Hamiltonian 
description with a superlattice potential, to show that the low energy theory contains an effective 
magnetic field. Numerical diagonalization of the Dirac equation confirms the presence of Landau 
levels. Possible consequences for transport are discussed. 
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Graphene is a honeycomb lattice of carbon atoms with 
a host of remarkable electronic properties [l|, 0] ■ Many 
of these are a result of the low energy states, which 
are governed by a massless Dirac equation. When un- 
doped the Fermi surface of the system is a set of dis- 
crete points, with resulting properties that compromise 
between metallic and insulating behavior. For example, 
an ideal graphene system can conduct diffusively even in 
the absence of impurities to scatter the electrons [!, 0] ■ 
The possibility of a universal conductivity even in the 
presence of impurities remains an active subject of discus- 
sion [H, H, [H~0] ■ The "relativistic" nature of the electronic 
spectrum introduces Klein paradox physics [l|, Q which 
may be evident in p — n junctions and related structures 
[9l-[l3|. Beyond this, graphene exhibits signatures of "ef- 
fective time-reversal symmetry breaking," effects most 
naturally described in terms of effective magnetic flux, 
due to a Berry's phase intrinsic to graphene [T^ - [l6j . or 
to strain in the lattice system [l7l - fl9| , even in the absence 
of a real applied magnetic field. 

In this paper we discuss a completely different way to 
induce an effective magnetic field in graphene, of suf- 
ficient uniformity that Landau levels in principle should 
be apparent in the electronic spectrum. This involves im- 
mersing the graphene system in a unidirectional periodic 
(i.e., superlattice) potential, for example by subjecting 
a sample with periodic ripples [2(| to a perpendicular 
electric field. Recently it has been shown [2l|, that 
such a potential can induce new Dirac points from the 
original one, whose positions arc determined by the am- 
plitude Vq of the periodic potential, and by its period. 
We demonstrate that by promoting Vq to a slowly vary- 
ing function of position, the effective theory for states 
near an induced Dirac point contains a gauge field, such 
that there are states in the spectrum corresponding to 
Landau levels. 

The possibility of generating such a description turns 
out to be intimately related to a chiral symmetry which is 
present in the Dirac equation even with a periodic poten- 
tial, as we show below. The symmetry requires that any 
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FIG. 1: (Color online) Energy bands for Dirac equation in a 
modulated superlattice potential. Inset: Ratio of first to sec- 
ond excited state energies near k y a/2n — 1.3, demonstrating 
the v2 ratio (dotted line) expected of Landau levels. 



state with energy e have a chiral partner at energy —e. 
These are degenerate for e — 0. Eigenstates near zero en- 
ergy can be approximately constructed from these chiral 
partners, and they have a number of properties guaran- 
teeing that the resulting effective theory has the form 
of a Dirac equation with a vector potential, representing 
a magnetic field with strength proportional to the rate 
at which the position of the induced Dirac point varies 
with position. We demonstrate the effect explicitly by 
numerically solving the Dirac equation with a potential 
Vo(x) cos Gqx. Results for the energy spectrum are il- 
lustrated in Fig. [TJ Within this spectrum we find level 
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spacings that are essentially identical to those of Landau 
levels (see inset of Fig. [I}, and wavefunctions contained 
in envelopes whose center positions - guiding centers - 
are fixed by a wavevector, and have the forms of har- 
monic oscillator states. We find that states that straddle 
regions in which dVo/dx changes sign are strongly dis- 
persive, acting much like edge states of graphene ribbons 
in a magnetic field [23| . The presence of such edge states 
should have very interesting consequences for transport 
in this system. 

Superlattice Potential and Chiral Symmetry - Low en- 
ergy states of graphene in a single valley may be modeled 
by a Dirac equation [23|. When an external potential is 
added, the eigenvalue equation to be solved has the form 
Hif) — sip, with H = vp(a x p x + u y p y ) + Vqv(x)I, where 
-0 is a two-component spinor whose entries are the wave- 
function amplitudes for the two sublattices, vf is the 
speed of electrons with momenta near the Dirac point, 



Pauli matrices, p x 



-ihd X y, Vq is the scale of 



an external potential with functional form v(x), and I is 
a unit 2x2 matrix. For Vq = 0, H is known to have a 
chiral symmetry, which can be expressed as {<J Zl H} = 0, 
with o~ z the third Pauli matrix. This anticommutator 
relation implies that for any eigenstate if) with energy e, 
there is a particle- hole partner o~ z ?p which is also an eigen- 
state of H, with energy — e. Moreover one may use this 
property to show that the state at zero energy is doubly 
degenerate. 

When there is a non-vanishing periodic superlattice 
potential satisfying v(x + a/2) = —v(x) with energy scale 
Vq, it was shown recently [U [22J that for sufficiently 
large V^a one obtains new Dirac cones in the spectrum. 
These appear as pairs which emerge from k = along 
the k y axis, one moving up and the other down as V$a 
increases, so that there are new zero energy states at 
k = (0,±A;*). Further increasing V$a generates yet more 
pairs of new Dirac cones. For simplicity, in what follows 
we will only consider the situation where a single new 
pair has emerged. 

The Hamiltonian with this type of periodic potential 
also has chiral symmetry, although the associated opera- 
tor is slightly more complicated. Defining a shift opera- 
tion Sip(x,y) — tp(x + a/2 7 y), it takes the form T = <r z S; 
one may easily demonstrate {H, T} = even in the pres- 
ence of the periodic potential. Because there is a chiral 
operator, one may show that a zero energy state (f>o at 
(0, k*) has a chiral partner Tcj>o which is also at zero en- 
ergy at the same wavevector. For what follows it is con- 
venient to define eigenstates of T, <j)± = [4>o ± T(j) ]/Af±, 
where Af± are normalization integrals. These states have 
a number of interesting properties which will be useful to 
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where purely real. Eq. Q] is easily shown 

by observing that <j)± are eigenstates of T with dif- 
ferent eigenvalues. Eq. [2] follows from {<J x , y ,o- z } — 
and < 0±|T 2 |</>± >= 1. Eqs. [3] and H follow 
from a specific form which </>o may take, 0o(r) — 
4>(x)[ismd(x),cos6(x)] 1i e lk y y , where ^ and 6 are real pe- 
riodic functions The quantities u x ^ y may be under- 
stood as the velocities of the electron states in the (x, y) 
directions of a Dirac cone, in units of vf, as we shall see 
below. 

Modulated Superlattice - We now consider the situation 
Vo — > Vb(x), where Vq varies slowly on the scale of the 
superlattice parameter a. The eigenfunctions generated 
from cj)Q have an extra spatial dependence from their im- 
plicit dependence on Vo- 4>o{x,y;Vo) — !• (f>o[x, y; Vq(x)]. 
In what follows it will be convenient to consider func- 
tions where the explicit spatial dependence and that 
inherited from Vq(x) are independent, so that we pro- 
mote the wavefunctions to functions of two variables: 
4>o[r';Vo(r)} — > <fio(r',r). Note that in this construction, 
if r is held constant then <fio(x' + a,y';r) = 4>o(x' ,y';r). 
We can then define matrix elements of the form 

px+a/2 />oo 

(^|A|0 2 ) r = / dx 1 dy'g(y-y')$ 1 {v'^A$ 2 {v\v) 

J x~ a/2 J -co 

where g{y— y') is a peaked function centered at and of 
width much greater than 2ir/k* [24| which integrates to 

1, and A is an operator. Identifying <fi± = 4>±(r'; Vq(x)), 
Eqs. [T]|4]remain true if we replace the standard matrix el- 
ements with those above. We now consider a generalized 
eigenvalue problem H(r', r)ip(r', r) = eV^f , r ); where 

H(r',r) = v F [-ia x (d x >+d x )-i<Ty(dy>+dy)]+Vo(x)v(x')2, 

and we have set H = 1. In the limit r' — > r, this becomes 
the eigenvalue equation one needs to solve for the physical 
Hamiltonian; the advantage of the generalized problem 
is that it retains the chiral symmetry {H,T} = 0, with 
the shift operation (S) contained in T acting only x' . To 
approximately solve this problem, we consider solutions 
of the form [24| 

^(r',r) = V+(r)<Mr',r) + ip-(r)$_(r' , r). 

Acting on this wavefunction with H, and then forming 
matrix elements (4>±\H — e\ip) r , we seek choices of the 
coefficients ip± for which these matrix elements vanish. 
After some algebra [25(, with the use of Eqs. [T]- SI one 
finds that they must obey the equations 

v x d x - ivydy + a+- (r) \ ( ip + 

-v x d x - ivydy + a_ + (r) J \ ip- 

Here V x>y {x) = VFU x , y [V (x)], and 

a±. T (r) = -iv F (<p±\v ■ V|(/> T ) r 
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It is important to note because the gradient operator in 
a± T acts on r and not r', a± ^ a^.. Thus the matrix 
appearing in Eq. [IT] is not Hermitian. However, one may 
show that lma± = — Ima± = v x yd x k* = v x A x is purely 
antisymmetric, while the real antisymmetric component 
of a, = Re(a± — a=p)/2, can be removed by a similarity 
transformation. Rewriting Eq. [17] in terms of ip± = 
exp[J dx a-^(x)vy(x)/v x (x)]ip±, we obtain a Hermitian 
eigenvalue equation Hip = eip, with 

H = (T y v x [-id x - A x ] + a x v y [-idy - Ay], (6) 

where A y — — Re(a± + a T )/2v y . Eq. [6] has the form of 
a Dirac Hamiltonian with a vector potential, albeit with 
the roles of the x and y momenta interchanged, and the 
effective velocities v XlV (x) being direction and position 
dependent. Several remarks are in order, (i) For Vo con- 
stant, v x , y lose their position dependence, and we expect 
them to give the velocities of the induced Dirac points for 
an unmodulated superlattice potential. Fig. [5] illustrates 
by comparison with velocities found numerically around 
an induced Dirac point of a graphene superlattice [U [22[ 
that this is indeed the case, (ii) Chiral symmetry plays 
a central role in giving Eq. [6] its Dirac form, through the 
properties of the <p± states expressed in Eqs. [T]|4] In par- 
ticular they guarantee that diagonal elements do not ap- 
pear in the Hamiltonian, and that the off-diagonal terms 
have a form that can naturally be grouped together into 
canonical momenta with vector potential contributions, 
(iii) The effective field represented by the vector potential 
is B e ff — d x A y — d y A x . In the limit of a slowly varying 
modulation (|a(9 x Vb)/Vo| <C 1), the first term is second 
order in derivatives and is negligibly small, so that the 
effective field is B eff = -d x k* y = -(dv k*[Vo})(d x V (x)). 
We see that the presence of an induced magnetic field is 
intimately connected with the fact that the position of 
the Dirac point in the Brillouin zone can be changed via 
an external parameter, (iv) With a gauge transforma- 
tion, H can be written in a form for which there is no 
explicit y dependence, so that k y is a good quantum num- 
ber, as we expect because the system is translationally 
invariant along the y direction. If we define a magnetic 
length 1/£q = d x k* then X = k y l^ becomes a guiding 
center coordinate around which the wavefunction is cen- 
tered. Approximating v XlV (x) —> v x , y (X), H becomes 
a standard Hamiltonian for Dirac electrons in a mag- 
netic field; the approximation is valid provided £q is suf- 
ficiently small. In this case one expects eigenenergies of 
the form e n (X) — sgn(n)yj2v x v y \n\B e f / with n an inte- 
ger. Although the values of e n {X) depend on k y through 
v x v y B e ff, neighboring energy levels n\, n-i at fixed k y 
should obey e ni / £« 2 = yWA*2- This is an unambiguous 
signature of Landau levels, which we next demonstrate 
numerically to be present in the energy spectrum of this 
system. 

Numerical Studies - Returning to the original Hamil- 
tonian, H — vp(a x p x + <JyP y ) + Vq(x)v(x)I, as a concrete 
example we take Vq(x) to have a sawtooth shape (illus- 
trated as a dashed line in Fig. [3]) , acting as an envelop for 
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FIG. 2: (Color online) Comparison of velocities v x , y , as com- 
puted from Eqs. [3j and [4] (solid lines), and from the nu- 
merically generated band structure in a superlattice potential 
(dots), around an induced Dirac point, as a function of Vo- 

a pure cosine potential v(x) = cosGox, with Go = 2n/a. 
The period of Vo is chosen to be Lq = 50a, and its mini- 
mum/maximum values are (l±0.1)[/o, with Uq = 0.12e_F 
and €f = h,VF ja. We expand H in plane wave states, for 
fixed values of k x (the Bloch wavevector) and k y , and 
diagonalize the resulting matrix. A typical example of 
a few of the eigenvalues as a function of k y , for k x = 0, 
are illustrated in Fig. [TJ Near k y — the Dirac point 
which would also be present for Vq — is apparent, but 
at larger \k y \ one finds bands that vary nearly linearly 
with k y , in the vicinity of the emergent Dirac points one 
would find if Vq(x) were replaced by its spatial average 
in the potential [2l|, |22j. These are the Landau levels 
generated by the potential modulation. The zero energy 
states one expects of Dirac Landau bands are very ap- 
parent. The inset of Fig. [T]shows the ratio of energies for 
the two lowest positive energy states; for Landau levels 
these should scale as yjnz/rii — \/2, which is indeed the 
case as illustrated in the inset. The levels shown are each 
very nearly doubly degenerate. This is because our cho- 
sen Vo(x) satisfies dV~o(x)/dx = — dVo(x + L /2)/dx, so 
that the effective magnetic field has regions of equal mag- 
nitude but opposite directions. We find that the spectra 
in the vicinity of the Landau levels are quantitatively 
nearly identical for any choice of k x . 

Fig. [3] illustrates wavefunction overlaps with cj)± for 
the lowest positive energy states for a single choice of k y . 
Their similarity to a lowest harmonic oscillator state for 
one of these, and to a first excited excited state for the 
other, is apparent. This is precisely what one expects for 
eigenstates of the Dirac equation in a uniform magnetic 
field [l[ . Note that the wavefunction has support in both 
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FIG. 3: (Color online) Wavefunctions for k y a/2n — 1.25, for 
a potential modulation Vo(x) of the form illustrated by the 
dashed line. 



the positive and negative effective field regions, but the 
roles of and <j>- which multiply the lower or higher 
harmonic oscillator states are reversed. This reflects the 
fact that Landau level raising and lowering operators in- 
terchange roles when the field sign is reversed. 

It is interesting to note that the degeneracy of the pos- 
itive and negative field states breaks down when k y is 
chosen so that the wavefunction is significantly differ- 
ent than zero at the cusps of Vq(x). In this situation 
the energy rises quickly as a function of k y , indicating 
a large current. These behave much as edge states in a 
quantum Hall system, and should have interesting con- 



sequences for transport. In particular we expect that 
o yy , the conductance along the y direction, should have 
steps as a function of Fermi energy Ep, occurring at en- 
ergies approximately reflecting the ^nB e f f behavior of 
the Landau levels. An observation of this effect would 
constitute evidence of the effective magnetic field present 
in this system. It is also interesting to note that left- and 
right-moving states are spatially separated, existing sep- 
arately at the upward and downward cusps of Vb (x) . We 
thus expect that backscattering due to disorder will be 
somewhat suppressed, although scattering into a time- 
reversed partner from the opposite valley will be present 
if disorder at short wavelengths is significant in the sys- 
tem. Studies of this and related transport phenomena 
are currently underway. 

In summary, we have demonstrated that an effective 
magnetic field can be generated in a graphene superlat- 
tice by a slow modulation of the superlattice potential 
amplitude. A description in terms of a gauge field can 
be produced from zero energy states at a Dirac point 
of the unmodulated superlattice which are chiral part- 
ners. These may be used as a basis to find states in the 
vicinity of the modulated potential. For a slowly varying 
amplitude, it was found that the effective field is pro- 
duced by the changing position of the Dirac point in real 
space. Numerical studies of the Dirac equation with this 
type of potential verified the presence of Landau levels in 
the spectrum, and suggest the possibility of interesting 
transport phenomena in this system. 
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Effective Magnetic Fields in Graphene Superlattices, Supplementary Information: 

Derivation of Equation 5 

In this supplement, we provide some details of how Eq. 5 is obtained from the generalized eigenvalue problem, 
H(r' ,r)t/j(r' ,r) = Eip(r',r), with 

H(r',r) = hv F [-ia x (d x , + d x ) - ia y (d y , + d y )] + V {x)v(x')X. (7) 
We consider approximate solutions of the form 

tp(r',r) = ■0+(r)0 + (r',r) + ip_(r)$_(r' , r), 
where (f>± = [<po ± T(/>o]/A/±, and Ho4>o = 0, with 

H (r',r) = hv F [-ia x d x > - icr y d y >} + V Q (x)v(x')l. (8) 

In this context, the operator T is defined as T<f>(r', r) = <j z S<fi(r' , r) = a z <fi(r' + ax/2, r). Since r is a fixed parameter, 
the potential in Hq is a periodic function of x' , and <j>o and <p± are zero energy Bloch states with respect to r'. As 
shown in Refs. 21 and 22, such states are present at k x = and k y = k*, with k* depending on Vq(x). Applying 

Eq. [5] to ip, multiplying by g(y — y')<j>±, and integrating r' over all y and a narrow interval of width a around x, one 
obtains two equations for the coefficients which may be expressed as 



(cf) + 1 5\<t>-) T ■ [—ihypVi^-] + [-ihv F (<p + \(j ■ V\c/>-)t\iP- + [—ihvp((j) + \a ■ v\(f>+) T \ip+ — eip+ 
(4>-\a\(j)+)r ■ [— ihvFVip+] + [—ifrvF(4>-\a ■ V|</>+) r ]V>+ + [— ihvp((f)-\a ■ V\<f>—) r ]ip- — eip—. 

In deriving this expression, we have used the fact that (4>±\<J XtV \(j)±) r = 0. 
We next demonstrate that (</>_|<? • V|</>_) r = 0. Writing 

' 1 <p-[x'; Vq(x)] ' 



(9) 



(10) 



we obtain explicitly 



x-\-a/2 poo 



ff-V|^_),= / dx' / g(y-y') 



x-a/2 



(11) 



The terms involving d x k* are proportional to J^^m dx'$- [r', Vq(:e)] ^a x 4>- [r', Vo(a;)], which vanishes identically, as 
in Eq. 2 of the text. What remains of this matrix element can be cast in the form 

rx+a/2 



\ff-V\<j>J) t = 



x a T d 



dx'4>-[x',y' = 0;F„(i)]V i S i ^[i',/ = 0;F 1 ( I )] 

x-a/2 
x+a/2 -i 

dx' - 

x-a/2 A/1 (a:) 

1 



<Mz';Vb(x)] - a z S4> [x';V (x)] 



Af-(x) 



<j)o[x'- : V (x)] - a z S<j> [x'; V (x)] 



(12) 



where in the last expression we have dropped the explicit y' in the arguments of the wave functions. Note the 
normalization factors A/"— obtain their spatial dependence from Vq(x). When the derivative operator acts on 1/A/"-, 
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the resulting contribution is again proportional to (<j>-\a x 



and vanishes. We now have 



|*-V|0_) r = 



x+a/2 



dx' 



x-a/2 



4>o [x'\ V (x)] - <j z S<j> [x'; V {x)} 



x a x d x 




dx'\ <p a x d x <f - [S(j>o^a x d x [S(fio\ 



(13) 



where we used {a x ,cr z } = and a z = 1. Finally, noting that S 1 is a Hermitian operator in x' and that S 2 = 1 since 
the wavefunctions are periodic in x', it easy to see that the four terms of the last expression mutually cancel. We 
thus have established (4>-\<J- V|^>_) r = 0. A very similar manipulation leads to (<j>+\cr • V|<^+) r = as well. Note that 
in reaching these results, the antisymmetry of a x and a z played a crucial role, so that the vanishing of this matrix 
element is in part a consequence of the fact that we can define a chiral operator for this system. That (<j>±\$ • V|</>±) r 
vanish is significant because had they not, the matrix equation for ip± would contain diagonal elements, effectively 
adding a scalar potential in addition to the vector potential we find generated by the modulation of Vq . 
Our equations now reduce to the form 



- ihv F (<f>+\cr x \(j>-) r d x ip- - ihv F (4> + \cr y \(f>-) r dy'il>- + -ihv F {4>+\5 ■ V\(f>-) r ip- = stp + 
-ihv F {(j)-\a x \(t) + ) T d x il) + - ihv F (<f>-\<j y \<f> + ) r dyip + + -iHv F {<j>-\<r ■ V\<f>+) r i>+ = eip- 

Using Eqs. 3 and 4 from the text, we can write four of the coefficients as position-dependent velocities: 

hv F ((t> + \(j x \(i)-) r = iv x (x), 
hv F ((f)^\a x \4> + ) r = -iv x (x), 

hV F ((j) + \ay\<j>-) r = Vy(x), 

hv F (<f>^\<7y\<t> + ) r = V y (x). 



(14) 



(15) 



Note again that the position dependence of these velocities enters through the x dependence of Vo(x). With the 
further definition 



a ±iT = -ihv F ((t>±\3 ■ V\4> T ) r = -i?iv F (<f)±\a x d x \<j> T )r, 



(16) 



where the second equality is a consequence of the fact that Vo depends on x but not on y, we arrive at Eq. 5 of the 
text: 



v x d x - ivydy + a + _(r) 

-v x d x - ivydy + a_+(r) 



